Using slimr to investigate the population genetics of the Sandy Inland Mouse in the periodic rainfall environment of the Simpson Desert
Russell Dinnage
30/04/2021
A look at the data
An ongoing project in the Simpson Desert of central Australia is producing genomic sequence data for small mammals and reptiles collected over more than 30 years. We have some preliminary data from this project for the Sandy Inland Mouse (Pseudomys hermannsburgensis) that we’d like to use to explore, through data analysis and simulation, and hopefully gain some insight into where this project might ultimately take us. First let’s load all the packages we use in the example, and then load the genetic data we will be working with. The data is stored in an .Rdata file containing a genlight object, which stores binary Single Nucleotide Polymorphism (SNP) data, defined in the adegenet package. We also have demographic data in the form of trap capture data (captures per 100 trap nights; see Dickman et al. 2014 and Greenville et al. 2016 for capture methods).
library(readr)
library(dplyr)
library(tidyr)
library(tibble)
library(mapview)
library(purrr)
library(conflicted)
library(future)
library(furrr)
library(lubridate)
library(ggplot2)
library(ggforce)
library(gganimate)
library(directlabels)
library(patchwork)
library(sf)
library(adegenet)
library(dartR)
library(here)
library(stringr)
library(slimr)
mapviewOptions(fgb = FALSE)
## because filter function exists in multiple packages we use
## the conflicted package to make sure we R uses the version in dplyr
conflict_prefer("filter", "dplyr")
conflict_prefer("initialize", "slimr")
## load data
gen <- read_rds("data/herm.rdata")
gen## /// GENLIGHT OBJECT /////////
##
## // 167 genotypes, 39,978 binary SNPs, size: 28.4 Mb
## 1028248 (15.4 %) missing data
##
## // Basic content
## @gen: list of 167 SNPbin
## @ploidy: ploidy of each individual (range: 2-2)
##
## // Optional content
## @ind.names: 167 individual labels
## @loc.names: 39978 locus labels
## @loc.all: 39978 alleles
## @position: integer storing positions of the SNPs
## @pop: population of each individual (group size range: 2-39)
## @other: a list containing: loc.metrics latlong ind.metrics loc.metrics.flags verbose historyabund <- read_csv("data/mammal_captures.csv")
abund## # A tibble: 1,408 x 8
## Year TripNo MonthYear SiteName SiteGrid `Notomys alexis` `Pseudomys herman~
## <dbl> <dbl> <chr> <chr> <chr> <dbl> <dbl>
## 1 1990 1 Mar.90 Main Camp MC1 3.70 0.463
## 2 1990 1 Mar.90 Main Camp MC2 1.67 1.11
## 3 1990 1 Mar.90 Main Camp MC4 1.85 0
## 4 1990 2 Jun.90 Main Camp MC1 2.38 0.397
## 5 1990 2 Jun.90 Main Camp MC2 1.19 1.19
## 6 1990 2 Jun.90 Main Camp MC4 1.59 1.19
## 7 1990 3 Aug.90 Main Camp MC1 0.397 0
## 8 1990 3 Aug.90 Main Camp MC2 0 2.38
## 9 1990 3 Aug.90 Main Camp MC4 0.794 0
## 10 1990 4 Sep.90 Main Camp MC1 1.39 4.17
## # ... with 1,398 more rows, and 1 more variable: Sminthopsis youngsoni <dbl>
So we have SNP sequences from 167 individuals taken over three years and several different sites. Before we take a look at the genetic data, let’s set the environmental context a bit. These small mammals live in the Simpson Desert, an arid ecosystem characterized by long periods of little rain, with occasional major rainfall events. These tend to occur every 6-10 years. After rainfall events, the desert turns green, it can be seen from space, and our analysis of satellite reflectance data confirms massive spikes in productivity after these rainfalls. Plants grow, flower and seed massively during these periods, and so food is abundant for small mammals. We can see the regular pulse of rainfall just by looking at the mammal abundances – Pseudomys hermannsburgensis in particular responds quickly to rainfall. Let’s look at our live trapping data on abundances to see this pattern.
## give more descriptive names to sites
sites <- tribble(~pop, ~SiteName,
"CS", "Carlo",
"FRN", "Field River North",
"FRS", "Field River South",
"KSE", "Kunnamuka Swamp East",
"MC", "Main Camp",
"SS", "South Site",
"WS", "Way Site"
)
abund_summ <- abund %>%
select(-`Notomys alexis`, -`Sminthopsis youngsoni`) %>%
left_join(sites) %>%
drop_na(pop) %>%
## convert month-year character column to proper dates using lubridate
separate(MonthYear, c("Month", "year"), "\\.") %>%
mutate(Month = sapply(Month, function(x) which(month.abb == x))) %>%
mutate(date = myd(paste(Month, year, sep = "-"), truncated = 1)) %>%
## summarise by site and date
group_by(date, pop) %>%
summarise(abund = mean(`Pseudomys hermannsburgensis`),
.groups = "drop") ## Joining, by = "SiteName"ggplot(abund_summ, aes(date, abund)) +
geom_path(aes(colour = pop)) +
geom_dl(aes(label = pop), method = "last.qp") +
aes(colour = pop) +
theme_minimal() +
theme(legend.position = "none")
Wow, those pulses are quite clear and quite extreme!
Okay, now let’s go have a look at the SNP data we have for these populations.
First, we will extract the metadata for the genetic samples and have a look at how they are distributed in the Simpson Desert.
gen_meta <- gen@other$ind.metrics
gen_meta <- gen_meta %>%
as_tibble() %>%
## fix some minor errors in genetic metadata with respect to site names:
mutate(pop = case_when(pop == "FR" ~ substr(as.character(SiteGrid), 1, 3),
pop == "KS" ~ "KSE",
TRUE ~ as.character(pop))) %>%
left_join(sites)## Joining, by = "pop"gen_meta## # A tibble: 167 x 18
## id pop id2 Date year month timem Species SiteGrid Trap Sex
## <chr> <chr> <int> <fct> <int> <int> <int> <fct> <fct> <fct> <fct>
## 1 P_herm~ KSE 121 16/05~ 2006 5 1 Pseudomys ~ KSEAS AS.32 "F"
## 2 P_herm~ FRN 193 17/09~ 2007 9 17 Pseudomys ~ FRN1 1.09 "F"
## 3 P_herm~ KSE 201 13/09~ 2007 9 17 Pseudomys ~ KSE1 1.28 "F"
## 4 P_herm~ SS 129 26/05~ 2006 5 1 Pseudomys ~ SS3 3.02 "M"
## 5 P_herm~ FRS 137 22/05~ 2006 5 1 Pseudomys ~ FRSU U.30 "M"
## 6 P_herm~ FRS 145 23/05~ 2006 5 1 Pseudomys ~ FRSU U.27 "M"
## 7 P_herm~ MC 153 28/05~ 2006 5 1 Pseudomys ~ MC2 2.31 "M"
## 8 P_herm~ CS 161 20/09~ 2007 9 17 Pseudomys ~ CS2G 2G.01 ""
## 9 P_herm~ CS 169 20/09~ 2007 9 17 Pseudomys ~ CSMTG MTG.~ ""
## 10 P_herm~ CS 177 21/09~ 2007 9 17 Pseudomys ~ CS2 2.03 ""
## # ... with 157 more rows, and 7 more variables: id_site <fct>, comments <fct>,
## # X <lgl>, lat <dbl>, lon <dbl>, datum <fct>, SiteName <chr>
Extract and plot site coordinates:
coords <- gen_meta %>%
group_by(pop) %>%
summarise(lon = mean(lon),
lat = mean(lat),
.groups = "drop") %>%
sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)
mapview(coords, label = coords$pop, map.types = "Esri.WorldImagery")Looks like our sites cluster into roughly 3 spatial groups. For simplicity we can merge these sites into 3 subpopulations, which will make it a bit simpler for our simulation later.
coords_3 <- gen_meta %>%
mutate(three_pop = case_when(pop %in% c("MC", "SS", "WS") ~ "BR",
pop %in% c("FRN", "FRS") ~ "BL",
pop %in% c("KSE", "CS") ~ "TR")) %>%
group_by(three_pop) %>%
summarise(lon = mean(lon),
lat = mean(lat),
.groups = "drop") %>%
ungroup() %>%
sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)
mapview(coords_3, label = coords_3$three_pop, map.types = "Esri.WorldImagery")As it happens these merged subpopulations are roughly equidistant, arranged in a triangle. If we assume the terrain between them is similarly easy to traverse we can simplify our model of the system even more and assume that all else being equal these subpopulations should have roughly similar migration rates between them. We called these subpopulations “BR”, “BL”, and “TR”, for Bottom Right, Bottom Left, and Top Right (creatively named I know).
Given this setup, let us have a look at the Fst values between these three subpopulations over the three years. We will use a custom function to calculate pairwise Fst for our subpopulations. We will do it this way so we can show some slimr functionality later, where we can exploit the similarity of R and SLiM code to run our Fst function from within SLiM, and get fast and live Fst estimates while a simulation is running. This means that the Fst calculated for our observed and our simulated data will be directly comparable when we do our downstream analysis.
sim.mutationFrequencies <- function(subpop) {
gen <- get("gen", parent.frame(2))
gl.alf(gen[pop = subpop])[ , 2]
}
isFinite <- function(x) is.finite(x)
calculateFST <- function(subpop1, subpop2) {
## Calculate the FST between two subpopulations
p1_p = sim.mutationFrequencies(subpop1);
p2_p = sim.mutationFrequencies(subpop2);
mean_p = (p1_p + p2_p) / 2.0;
H_t = 2.0 * mean_p * (1.0 - mean_p);
H_s = p1_p * (1.0 - p1_p) + p2_p * (1.0 - p2_p);
fst = 1.0 - H_s/H_t;
fst = fst[isFinite(fst)]; ## exclude muts where mean_p is 0.0 or 1.0
return(mean(fst));
}
pairFST <- function(gen) {
pops <- as.character(unique(pop(gen)))
pop_pairs <- combn(pops, 2)
pair_fst <- apply(pop_pairs, 2, function(x) calculateFST(x[1], x[2]))
fst_mat <- matrix(nrow = length(pops), ncol = length(pops))
rownames(fst_mat) <- colnames(fst_mat) <- pops
fst_mat[t(pop_pairs)] <- pair_fst
t(fst_mat)
}We used several custom functions above and named them so that they would match available functions in SLiM. That way, when we run our simulations using slimr, SLiM will use its optimized versions of these functions internally. We also create a function to run the pairwise Fst function on all possible subpopulation pairs, which we will now run on our observed SNP data. We used a function in the dartR package to calculate allele frequencies to use in our Fst calculation. In SLiM, this is a builtin function, which we will take advantage of later. Note that the use of get and parent.frame is a trick to get R to use the right gen object without having to specify it as an argument, which will let us reuse the functions later in SLiM, which cannot access R objects while it is running.
gen_meta <- gen_meta %>%
mutate(three_pop = case_when(pop %in% c("MC", "SS", "WS") ~ "BR",
pop %in% c("FRN", "FRS") ~ "BL",
pop %in% c("KSE", "CS") ~ "TR"))
gen@other$ind.metrics <- gen_meta
pop(gen) <- gen_meta$three_pop
## remove population we have genetic data for but no abundance data
gen <- gen[gen@other$ind.metrics$pop != "WS"]
fst <- pairFST(gen)
fst## TR BL BR
## TR NA NA NA
## BL 0.01437133 NA NA
## BR 0.01358429 0.0156385 NAThat is the overall Fst, but we really want to know how that is changing over our three sampling time points. We have three years of data, which happen to span a major rainfall event, from 2006, 2007, and 2008, where there was a major rainfall event in 2007. Let’s look at the pairwise Fst over the three years.
fst_by_year <- map(c(2006, 2007, 2008),
~pairFST(gen[gen@other$ind.metrics$year == .x, ]) %>%
as.dist() %>%
as.matrix())
names(fst_by_year) <- c("2006", "2007", "2008")
fst_by_year## $`2006`
## TR BR BL
## TR 0.00000000 0.07423650 0.07492284
## BR 0.07423650 0.00000000 0.04449832
## BL 0.07492284 0.04449832 0.00000000
##
## $`2007`
## BL TR BR
## BL 0.00000000 0.02405952 0.03053511
## TR 0.02405952 0.00000000 0.02336933
## BR 0.03053511 0.02336933 0.00000000
##
## $`2008`
## BL BR TR
## BL 0.00000000 0.06184445 0.06037648
## BR 0.06184445 0.00000000 0.04472136
## TR 0.06037648 0.04472136 0.00000000We can visualize that over the three years:
fst_df <- imap_dfr(fst_by_year,
~combn(c("BR", "BL", "TR"), 2) %>%
t() %>%
as.data.frame() %>%
rename(pop1 = V1, pop2 = V2) %>%
mutate(fst = .x[cbind(pop1, pop2)],
year = .y,
pop_combo = paste(pop1, pop2, sep = " to ")))
ggplot(fst_df, aes(year, fst)) +
geom_path(aes(colour = pop_combo, group = pop_combo)) +
theme_minimal()
So we see that during the rainfall year a precipitous drop in Fst occurs, which has only partially recovered the following year. Now, can we think of a reason why this pattern might occur? We have a hypothesis! First, we’ll lay it out, then we will interrogate its logic using the simulation tools provided by slimr. Earlier on we talked about how the desert “turns green” after a big rainfall event. This implies it reverts back to dry or “red” sands afterwards. Generally this is true, however, some patches remain relatively green during the dry periods. These are sometimes know as mesic ‘refugia’. We hypothesized that during dry periods, our mouse population recedes into these separate refugia across the landscape, where they become relatively reproductively isolated. This leads to an inexorable increase in Fst due to a genetic drift without gene flow. Then, when the rain comes, the whole landscape becomes wet, the mice come out to play, spreading throughout, reproducing along the way. The subpopulations of the refugia intermix, leading to an erasure of the genomic differentiation that had been in progress during the dry years. But can this mechanism lead to the level of observed change in such a short period? Are there other possible mechanisms? Let’s explore these ideas with simulation!
A simple simulation
To begin setup a simple simulation that captures some of the features of the system we are studying. We will simulate three subpopulations, as we have in our data, and allow migration between them at equal rates. We will include evolutionary processes such as mutation, selection, sexual reproduction, and recombination as additional processes happening in our simulation. Our initial question is simply, can we observe Fst values that change through time in a manner consistent with our observed changes using combinations of just these processes? Thanks to slimr, we can continue and write out the logic of our proposed simulation directly in R, using syntax very similar to SLiM (>3.0). We can additionally make the simulation code easy to update with different parameter combinations to try out.
We will use the following parameters in our simulation:
| Parameter | Description |
|---|---|
| mut_rate | mutation rate of simulated population |
| genome_size | size of simulated genome |
| selection_strength | mean strength of selection specified as a standard deviation of selection coeffs |
| migration_rates | rate of migration amongst three pops when abundance is high, must be between 0 and 1 |
| abund_threshold | abundance (before scaling) above which migration between populations is “turned” on |
| recomb_rate | recombination rate, a nuisance parameter |
| popsize_scaling | multiply observed abundances by this value to get total subpop size |
Now for the code:
default_genome_size <- 300000
pop_sim <- slim_script(
slim_block(initialize(), {
setSeed(12345);
initializeMutationRate(slimr_template("mut_rate", 1e-6));
initializeMutationType("m1", 0.5, "n", 0, slimr_template("selection_strength", 0.1));
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, slimr_template("genome_size", !!default_genome_size) - 1);
initializeRecombinationRate(slimr_template("recomb_rate", 1e-8));
initializeSex("A");
defineConstant("abund", slimr_inline(pop_abunds, delay = TRUE));
defineConstant("sample_these", slimr_inline(sample_these, delay = TRUE));
}),
slim_block(1, {
init_pop = slimr_inline(init_popsize, delay = TRUE)
## set populations to initial size
sim.addSubpop("p1", asInteger(init_pop[0]));
sim.addSubpop("p2", asInteger(init_pop[1]));
sim.addSubpop("p3", asInteger(init_pop[2]));
}),
slim_block(1, late(), {
## get starting population from a file which we will fill-in later
sim.readFromPopulationFile(slimr_inline(starting_pop, delay = TRUE));
## migration on or off flags for pops 1-3 (using tag)
p1.tag = 0;
p2.tag = 0;
p3.tag = 0;
}),
slim_block(1, 1000, late(), {
## update generation number
gen = sim.generation %% 50
if(gen == 0) {
gen = 50
}
## set population size to observed levels
p1.setSubpopulationSize(asInteger(ceil(abund[0, gen - 1] * slimr_template("popsize_scaling", 100))));
p2.setSubpopulationSize(asInteger(ceil(abund[1, gen - 1] * ..popsize_scaling..)));
p3.setSubpopulationSize(asInteger(ceil(abund[2, gen - 1] * ..popsize_scaling..)));
## increase migration when above abundance threshold
if(p1.tag == 0 & abund[0, gen - 1] > slimr_template("abund_threshold", 5)) {
p2.setMigrationRates(p1, slimr_template("migration_rate", 0))
p3.setMigrationRates(p1, ..migration_rate..)
p1.tag = 1;
}
if(p1.tag == 1 & abund[0, gen - 1] <= ..abund_threshold..) {
p2.setMigrationRates(p1, 0)
p3.setMigrationRates(p1, 0)
p1.tag = 0;
}
if(p2.tag == 0 & abund[1, gen - 1] > ..abund_threshold..) {
p1.setMigrationRates(p2, ..migration_rate..)
p3.setMigrationRates(p2, ..migration_rate..)
p2.tag = 1;
}
if(p2.tag == 1 & abund[1, gen - 1] <= ..abund_threshold..) {
p1.setMigrationRates(p2, 0)
p3.setMigrationRates(p2, 0)
p2.tag = 0;
}
if(p3.tag == 0 & abund[2, gen - 1] > ..abund_threshold..) {
p1.setMigrationRates(p3, ..migration_rate..)
p2.setMigrationRates(p3, ..migration_rate..)
p3.tag = 1;
}
if(p3.tag == 1 & abund[2, gen - 1] <= ..abund_threshold..) {
p1.setMigrationRates(p3, 0)
p2.setMigrationRates(p3, 0)
p3.tag = 0;
}
}),
slim_block(1000, late(), {
slimr_output_full()
})
)
pop_sim## <slimr_script[5]>
## block_init:initialize() {
## setSeed(12345);
## initializeMutationRate(..mut_rate..);
## initializeMutationType("m1", 0.5, "n", 0, ..selection_strength..);
## initializeGenomicElementType("g1", m1, 1);
## initializeGenomicElement(g1, 0, ..genome_size.. - 1);
## initializeRecombinationRate(..recomb_rate..);
## initializeSex("A");
## defineConstant("abund", slimr_inline(pop_abunds));
## defineConstant("sample_these", slimr_inline(sample_these));
## }
##
## block_2:1 early() {
## init_pop = slimr_inline(init_popsize);
## sim.addSubpop("p1", asInteger(init_pop[0]));
## sim.addSubpop("p2", asInteger(init_pop[1]));
## sim.addSubpop("p3", asInteger(init_pop[2]));
## }
##
## block_3:1 late() {
## sim.readFromPopulationFile(slimr_inline(starting_pop));
## p1.tag = 0;
## p2.tag = 0;
## p3.tag = 0;
## }
##
## block_4:1:1000 late() {
## gen = sim.generation%%50;
## if (gen == 0) {
## gen = 50;
## }
## p1.setSubpopulationSize(asInteger(ceil(abund[0, gen - 1] * ..popsize_scaling..)));
## p2.setSubpopulationSize(asInteger(ceil(abund[1, gen - 1] * ..popsize_scaling..)));
## p3.setSubpopulationSize(asInteger(ceil(abund[2, gen - 1] * ..popsize_scaling..)));
## if (p1.tag == 0 & abund[0, gen - 1] > ..abund_threshold..) {
## p2.setMigrationRates(p1, ..migration_rate..);
## p3.setMigrationRates(p1, ..migration_rate..);
## p1.tag = 1;
## }
## if (p1.tag == 1 & abund[0, gen - 1] <= ..abund_threshold..) {
## p2.setMigrationRates(p1, 0);
## p3.setMigrationRates(p1, 0);
## p1.tag = 0;
## }
## if (p2.tag == 0 & abund[1, gen - 1] > ..abund_threshold..) {
## p1.setMigrationRates(p2, ..migration_rate..);
## p3.setMigrationRates(p2, ..migration_rate..);
## p2.tag = 1;
## }
## if (p2.tag == 1 & abund[1, gen - 1] <= ..abund_threshold..) {
## p1.setMigrationRates(p2, 0);
## p3.setMigrationRates(p2, 0);
## p2.tag = 0;
## }
## if (p3.tag == 0 & abund[2, gen - 1] > ..abund_threshold..) {
## p1.setMigrationRates(p3, ..migration_rate..);
## p2.setMigrationRates(p3, ..migration_rate..);
## p3.tag = 1;
## }
## if (p3.tag == 1 & abund[2, gen - 1] <= ..abund_threshold..) {
## p1.setMigrationRates(p3, 0);
## p2.setMigrationRates(p3, 0);
## p3.tag = 0;
## }
## }
##
## block_5:1000 late() {
## {sim.outputFull() -> full_output}
## }
## This slimr_script has templating in block(s) block_init
## and block_4 for variables mut_rate and
## selection_strength and genome_size and recomb_rate and
## popsize_scaling and abund_threshold and migration_rate.
The above script uses a number of different slimr feature that we will explain here. To anyone who has used SLiM before, the above script probably looks pretty familiar. For those of you not familiar with SLiM, we will go over each part of the above script one by one to gain an understanding of what the script is doing. Otherwise, if you’d like to know more about SLiM, I would highly recommend reading the SLiM manual, which is very detailed and packed full of examples.
We start by setting a variable that we will use in the slimr_script and later on in the R script as well.
We then use the slimr_script() function to create a new slimr_script object, which will contain all the information needed to run it in SLiM later. slimr_script takes a list of slimr_block objects as arguments, which are created with the slimr_block() function. slimr_block objects contain blocks of SLiM code. Generally, there are two main types of code blocks in SLiM: initialize blocks and ‘event’ blocks. Out first block above is an initialize block, which sets up the SLiM simulation. All slimr_script() calls must contain at least one initialize block. Most of the action with slimr happens in this block for this script. In it, we set up the simulation with our parameters, using two ‘slimr verbs’. ‘slimr verbs’ are R functions that can be inserted into SLiM code and which modify the code for various useful purposes. Here we are using slimr_template(), which lets you created ‘templated’ variables in a script – these are placemarkers that can be filled in with desired values later. We also use slimr_inline(), which let’s you ‘inline’ R objects directly into a script. Inlining means the object and its values are converted into a text representation and passed to SLiM by embedding them directly into the script. This all happens behind the scenes, so you don’t need to worry about it if you are not interested. We have used the special argument delay = TRUE to delay the evaluation of slimr_inline() until later. This is because we haven’t created the R objects slimr_inline() references yet.
The initialize block is followed by 4 more slimr_block calls, each of which specifies SLiM ‘events’. “block_2” sets up 3 subpopulations in generation 1, and “block_3” sets up starting conditions for these subpopulation by reading from a file (which we create later), also in generation 1. “block_4” contains most of the logic of the simulation, running in every generation from 1 to 1000. Finally, “block_5” outputs the state of the simulation at the last generation, generation 1000. We assigned the slimr_script object to a variable called pop_sim, and printed this object, obtaining a pretty output of our script, including special syntax showing our templated variable, which look like this: ..variable_name.. So, in order to run this script in SLiM, we first need to fill in values for our templated variables, and also make sure we the R objects referred to in slimr_inline() calls exist in our R session. We can then “render” our script into a runnable form.
so, how do we fill in the templated variables with values and fill-in our slimr_inline calls with actual R objects? For that, we use the slim_script_render function. The first thing we need to do before we use it is to create the R objects that slimr_inline refers to, otherwise we will get an error about non-existent objects. So what do we need? We need a matrix of abundances for our three subpopulations (for slimr_inline(pop_abunds)), a matrix of starting population abundances (for slimr_inline(init_popsize)) and a file name pointing to a SLiM population data file containing initial population conditions (for slimr_inline(starting_pop). We will later create this file from our genlight object using slimr. But first, let’s get our objects in order, choose some parameter values and render ourselves a slimr_script!
For our population abundances over time, we will use the live-trapping data, and we will create a simplified cycle, taking only part of the sequence, and then repeating it over and over throughout the years. We will sample data between 1995 and 2009, interpolate that to 50 time points, then loop over it inside our slimr_script, as we discussed above.
pop_abunds <- abund_summ %>%
filter(date < "2009-10-01" & date > "1995-03-01") %>%
mutate(three_pop = case_when(pop %in% c("MC", "SS", "WS") ~ "BR",
pop %in% c("FRN", "FRS") ~ "BL",
pop %in% c("KSE", "CS") ~ "TR")) %>%
drop_na(three_pop) %>%
group_by(date, three_pop) %>%
summarise(abund = mean(abund, na.rm = TRUE),
.groups = "drop")
pop_abunds <- pop_abunds %>%
pivot_wider(names_from = date, values_from = abund) %>%
as.matrix()
pop_abunds <- rbind(approx(pop_abunds[1, ], n = 50)$y,
approx(pop_abunds[2, ], n = 50)$y,
approx(pop_abunds[3, ], n = 50)$y)## Warning in xy.coords(x, y, setLab = FALSE): NAs introduced by coercion
## Warning in xy.coords(x, y, setLab = FALSE): NAs introduced by coercion
## Warning in xy.coords(x, y, setLab = FALSE): NAs introduced by coercion## replace exact zeroes
pop_abunds[pop_abunds == 0] <- 0.02
## set sample times corresponding to our genetic data (roughly)
sample_times <- c("2006" = 40, "2007" = 45, "2008" = 49)
## plot our abundance sequence
plot(pop_abunds[2, ], type = "l", col = "blue")
lines(pop_abunds[1, ], col = "red")
lines(pop_abunds[3, ], col = "green")
abline(v = sample_times)
So, that is out pop_abunds. For initial population states we are going to base it on data from 2008, which is the third sampling date, and which corresponds to what we think is a panmictic population, which is just after, but not during the big rainfall event. We will get the initial population sizes from here, as well as our starting population data. Our initial population sizes need to be based on this (rather than our abundance data we just generated above), because the number of individuals at the start of our simulation must match the number in the starting population data. The population will then be immediately adjusted to our desired population size by essentially drawing offspring from our smaller starting population pool. This should give us populations with SNP frequencies resembling our actual data at the start of the simulation.
## extract 2008 data
gen_2008 <- gen[gen@other$ind.metrics$year == 2008, ]
## count number of individuals in genetic sample per subpopulation
init_popsize <- c(table(pop(gen_2008)))
## set filename to be used for starting pop data (using slim_file to make sure SLiM can find it)
starting_pop = here("sims/starting_pop.txt") %>%
slim_file()
## setup generations to sample
## we will just sample the three years corresponding to the data, but do it during the last 6 cycles, as "technical" replicates
sample_these <- purrr::map(c(14:19),
~50*.x + sample_times) %>%
purrr::reduce(union)Okay, now we are ready to use slim_script_render to generate a script we can run. The first thing we can try is to render the script without providing a template. Since we have provided defaults for all of our templated variables, this should generate a script with the default values (our ‘default’ script).
script_def <- slim_script_render(pop_sim)## Warning: Warning: A templated variable was not specified in the template and has been replaced by its default value.script_def## <slimr_script[5]>
## block_init:initialize() {
## setSeed(12345);
## initializeMutationRate(1e-06);
## initializeMutationType("m1", 0.5, "n", 0, 0.1);
## initializeGenomicElementType("g1", m1, 1);
## initializeGenomicElement(g1, 0, 3e+05 - 1);
## initializeRecombinationRate(1e-08);
## initializeSex("A");
## defineConstant("abund", {. = matrix(c(0.926, 0.278, 0.463, 0.557, 0.153, ...), nrow = 3, ncol = 50)});
## defineConstant("sample_these", {. = c(740, 745, 749, 790, 795, 799, 840, 845, 849, 890, ...)});
## }
##
## block_2:1 early() {
## init_pop = {. = c(7, 13, 13)}
## sim.addSubpop("p1", asInteger(init_pop[0]));
## sim.addSubpop("p2", asInteger(init_pop[1]));
## sim.addSubpop("p3", asInteger(init_pop[2]));
## }
##
## block_3:1 late() {
## sim.readFromPopulationFile({. = "/mnt/f/Projects/slimr_resources/sims/starting_pop.txt"});
## p1.tag = 0;
## p2.tag = 0;
## p3.tag = 0;
## }
##
## block_4:1:1000 late() {
## gen = sim.generation%%50;
## if (gen == 0) {
## gen = 50;
## }
## p1.setSubpopulationSize(asInteger(ceil(abund[0, gen - 1] * 100)));
## p2.setSubpopulationSize(asInteger(ceil(abund[1, gen - 1] * 100)));
## p3.setSubpopulationSize(asInteger(ceil(abund[2, gen - 1] * 100)));
## if (p1.tag == 0 & abund[0, gen - 1] > 5) {
## p2.setMigrationRates(p1, 0);
## p3.setMigrationRates(p1, 0);
## p1.tag = 1;
## }
## if (p1.tag == 1 & abund[0, gen - 1] <= 5) {
## p2.setMigrationRates(p1, 0);
## p3.setMigrationRates(p1, 0);
## p1.tag = 0;
## }
## if (p2.tag == 0 & abund[1, gen - 1] > 5) {
## p1.setMigrationRates(p2, 0);
## p3.setMigrationRates(p2, 0);
## p2.tag = 1;
## }
## if (p2.tag == 1 & abund[1, gen - 1] <= 5) {
## p1.setMigrationRates(p2, 0);
## p3.setMigrationRates(p2, 0);
## p2.tag = 0;
## }
## if (p3.tag == 0 & abund[2, gen - 1] > 5) {
## p1.setMigrationRates(p3, 0);
## p2.setMigrationRates(p3, 0);
## p3.tag = 1;
## }
## if (p3.tag == 1 & abund[2, gen - 1] <= 5) {
## p1.setMigrationRates(p3, 0);
## p2.setMigrationRates(p3, 0);
## p3.tag = 0;
## }
## }
##
## block_5:1000 late() {
## {sim.outputFull() -> full_output}
## }
Okay, so now we can see that we have an indication that our R objects will be inserted into the SLiM script. We can try and run this in SLiM now, but we will get an error:
test <- slim_run(script_def)##
##
## Simulation finished with exit status: 1## Warning:
## Failed! Error:## ERROR (SLiMSim::InitializePopulationFromFile): initialization file does not exist or is empty.
##
## Error on script line 26, character 8:
##
## sim.readFromPopulationFile("/mnt/f/Projects/slimr_resources/sims/starting_pop.txt");
## ^^^^^^^^^^^^^^^^^^^^^^This is because the starting_pop.txt doesn’t exist yet, and SLiM has thrown an error telling us this (which slimr passes from SLiM to us through the R console). Let’s create the initial population file using slim_make_pop_init. We need to pass this function either a genlight or a SNP matrix. Since our genlight object has some missing values, and SLiM cannot handle missing values, we will first convert to a SNP matrix, and then fill-in missing values by interpolation. Essentially we just fill-in missing values randomly with a draw from a binomial distribution.
snp_mat <- as.matrix(gen_2008)
glimpse(snp_mat)## int [1:33, 1:39978] 0 0 0 0 0 0 0 0 0 0 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:33] "P_herm_281" "P_herm_289" "P_herm_297" "P_herm_265" ...
## ..$ : chr [1:39978] "15067413-48-G/A" "15067414-5-C/T" "15067416-63-G/A" "15067419-10-G/C" ...## replace NAs
nas <- apply(snp_mat, 2, function(x) !any(!is.finite(x)))
snp_mat[ , !nas] <- apply(snp_mat[ , !nas], 2, function(x) {
tab <- table(x[is.finite(x)]);
x[!is.finite(x)] <- as.integer(sample(names(tab), sum(!is.finite(x)), replace = TRUE, prob = tab / sum(tab)));
x
})
glimpse(snp_mat)## int [1:33, 1:39978] 0 0 0 0 0 0 0 0 0 0 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:33] "P_herm_281" "P_herm_289" "P_herm_297" "P_herm_265" ...
## ..$ : chr [1:39978] "15067413-48-G/A" "15067414-5-C/T" "15067416-63-G/A" "15067419-10-G/C" ...Then, we make our starting population file!
sexes <- as.character(gen_2008@other$ind.metrics$Sex)
## sex ratio of sample is skewed to female, so replace actually sexes with extra males
## otherwise simulation crashes because it is unable to sample enough females? Need
## to look more into this bug.
sexes[pop(gen_2008) == "BL"] <- c("M", "F")## Warning in sexes[pop(gen_2008) == "BL"] <- c("M", "F"): number of items to
## replace is not a multiple of replacement lengthsexes[pop(gen_2008) == "BR"] <- c("M", "F")## Warning in sexes[pop(gen_2008) == "BR"] <- c("M", "F"): number of items to
## replace is not a multiple of replacement lengthsexes[pop(gen_2008) == "TR"] <- c("M", "F")## Warning in sexes[pop(gen_2008) == "TR"] <- c("M", "F"): number of items to
## replace is not a multiple of replacement length## now make initial population file
slim_make_pop_input(snp_mat, here("sims/starting_pop.txt"), ## filename
sim_gen = 1, ## set generation to first generation
ind_pops = gen_2008@other$ind.metrics$three_pop, ## use subpops
ind_sex = sexes, ## set sexes
## random positions:
mut_pos = sample.int(default_genome_size - 1, nLoc(gen_2008)),
version = 3)## [1] "F:/Projects/slimr_resources/sims/starting_pop.txt"## look at the first 50 line to see if it worked:
read_lines(here("sims/starting_pop.txt")) %>%
head(50)## [1] "#OUT 1 A" "Version: 3"
## [3] "Populations:" "p1 7 S 0.42857142857142855"
## [5] "p2 13 S 0.46153846153846156" "p3 13 S 0.46153846153846156"
## [7] "Mutations:" "0 0 m1 28309 0 0.5 p1 1 2"
## [9] "1 1 m1 43642 0 0.5 p1 1 8" "2 2 m1 47100 0 0.5 p1 1 4"
## [11] "3 3 m1 16002 0 0.5 p1 1 19" "4 4 m1 9943 0 0.5 p1 1 8"
## [13] "5 5 m1 298798 0 0.5 p1 1 4" "6 6 m1 18014 0 0.5 p1 1 6"
## [15] "7 7 m1 16143 0 0.5 p1 1 5" "8 8 m1 297437 0 0.5 p1 1 2"
## [17] "9 9 m1 46999 0 0.5 p1 1 2" "10 10 m1 140145 0 0.5 p1 1 8"
## [19] "11 11 m1 71064 0 0.5 p1 1 10" "12 12 m1 220419 0 0.5 p1 1 26"
## [21] "13 13 m1 1016 0 0.5 p1 1 2" "14 14 m1 290905 0 0.5 p1 1 1"
## [23] "15 15 m1 227173 0 0.5 p1 1 66" "16 16 m1 91797 0 0.5 p1 1 3"
## [25] "17 17 m1 88114 0 0.5 p1 1 66" "18 18 m1 175169 0 0.5 p1 1 2"
## [27] "19 19 m1 188909 0 0.5 p1 1 62" "20 20 m1 71616 0 0.5 p1 1 12"
## [29] "21 21 m1 29804 0 0.5 p1 1 17" "22 22 m1 201396 0 0.5 p1 1 64"
## [31] "23 23 m1 256687 0 0.5 p1 1 6" "24 24 m1 73500 0 0.5 p1 1 2"
## [33] "25 25 m1 178485 0 0.5 p1 1 3" "26 26 m1 166516 0 0.5 p1 1 14"
## [35] "27 27 m1 283175 0 0.5 p1 1 18" "28 28 m1 51058 0 0.5 p1 1 1"
## [37] "29 29 m1 290828 0 0.5 p1 1 1" "30 30 m1 94668 0 0.5 p1 1 1"
## [39] "31 31 m1 189202 0 0.5 p1 1 6" "32 32 m1 78166 0 0.5 p1 1 3"
## [41] "33 33 m1 214841 0 0.5 p1 1 65" "34 34 m1 80153 0 0.5 p1 1 1"
## [43] "35 35 m1 89149 0 0.5 p1 1 1" "36 36 m1 227494 0 0.5 p1 1 20"
## [45] "37 37 m1 141882 0 0.5 p1 1 12" "38 38 m1 426 0 0.5 p1 1 18"
## [47] "39 39 m1 13410 0 0.5 p1 1 3" "40 40 m1 63745 0 0.5 p1 1 63"
## [49] "41 41 m1 57925 0 0.5 p1 1 11" "42 42 m1 167901 0 0.5 p1 1 9"And now we can try running our simulation:
res <- slim_run(script_def, throw_error = TRUE)##
##
## Simulation finished with exit status: 0##
## Success!Okay, so now we have a working simulation! We’ve currently just output the full population information at the end of the simulation, which doesn’t give us too much information about our question, but just to convince ourselves the simulation produced genetic information, we can extract our output as a genlight object and plot it. Note: The conversion can take awhile.
pop_res <- slimr::slim_extract_genlight(res, "full_output")
## reorder by subpop
pop_res <- pop_res[order(pop(pop_res)), ]
plot(pop_res)
So this shows that our three populations are characterized by mostly non-overlapping fixed mutations. This makes sense since our default parameters specified some selection, and no migration. Mutation with positive selection coefficients would have sweeped to fixation, and the lack of migration between subpopulation means that each subpopulation is likely to have had a different set of mutations go to fixation, and no way for these to then spread to other subpopulations.
Okay, to get an idea of what sort of dynamics of Fst we will get with different versions of this simulation, let’s create a version where we calculate Fst during the simulation, and output it for tracking. Then we can make some time-series plots. When it comes to comparing simulation outputs to our actual data we will make yet another version where we will only output just before, during, and just after the population explosions corresponding to rainfall. To add a Fst calculation to our simulation we will use the slim_function function, which allows you to specify a SLiM function to use in your SLiM code. This Fst function is based on the one described in the SLiM manual, and can also be found online here: https://github.com/MesserLab/SLiM-Extras/blob/master/functions/calcFST.txt. We have already used the function as an R function when we calculated Fst for our observed data. Now we can convert it to a SLiM function. slim_function requires us to specify the function arguments and the function names as character strings, and then the body of the function as slimr code. If inserted into a slim_script call, the function will then be accessible inside the rest of the code in the slim_script call. We can use the previously specified R function code body by simply extracting it using rlang::fn_body(). We also need to make sure it is evaluated before insertion, which we assure using the evaluation forcing operator !!.
Note: The latest version of SLiM (v3.5), which was released after the initial version of this vignette was written, now has a built-in function called calcFST, which means we technically no longer need to use a custom function for this functionality. However, we have kept this section of the vignette to show how custom functions can be specified and used in slimr.
slim_script(
slim_function("o<Subpopulation>$ subpop1", "o<Subpopulation>$ subpop2",
name = "calculateFST",
return_type = "f$", body = calculateFST),
slim_block(initialize(), {
setSeed(12345);
initializeMutationRate(slimr_template("mut_rate", 1e-6));
initializeMutationType("m1", 0.5, "n", 0, slimr_template("selection_strength", 0.1));
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, slimr_template("genome_size", !!default_genome_size) - 1);
initializeRecombinationRate(slimr_template("recomb_rate", 1e-8));
initializeSex("A");
defineConstant("abund", slimr_inline(pop_abunds, delay = TRUE));
defineConstant("sample_these", slimr_inline(sample_these, delay = TRUE));
}),
slim_block(1, {
init_pop = slimr_inline(init_popsize, delay = TRUE)
## set populations to initial size
sim.addSubpop("p1", asInteger(init_pop[0]));
sim.addSubpop("p2", asInteger(init_pop[1]));
sim.addSubpop("p3", asInteger(init_pop[2]));
}),
slim_block(1, late(), {
## get starting population from a file which we will fill-in later
sim.readFromPopulationFile(slimr_inline(starting_pop, delay = TRUE));
## migration on or off flags for pops 1-3 (using tag)
p1.tag = 0;
p2.tag = 0;
p3.tag = 0;
}),
slim_block(1, 1000, late(), {
## update generation number
gen = sim.generation %% 50
if(gen == 0) {
gen = 50
}
## set population size to observed levels
p1.setSubpopulationSize(asInteger(ceil(abund[0, gen - 1] * slimr_template("popsize_scaling", 100))));
p2.setSubpopulationSize(asInteger(ceil(abund[1, gen - 1] * ..popsize_scaling..)));
p3.setSubpopulationSize(asInteger(ceil(abund[2, gen - 1] * ..popsize_scaling..)));
## increase migration when above abundance threshold
if(p1.tag == 0 & abund[0, gen - 1] > slimr_template("abund_threshold", 5)) {
p2.setMigrationRates(p1, slimr_template("migration_rate", 0))
p3.setMigrationRates(p1, ..migration_rate..)
p1.tag = 1;
}
if(p1.tag == 1 & abund[0, gen - 1] <= ..abund_threshold..) {
p2.setMigrationRates(p1, 0)
p3.setMigrationRates(p1, 0)
p1.tag = 0;
}
if(p2.tag == 0 & abund[1, gen - 1] > ..abund_threshold..) {
p1.setMigrationRates(p2, ..migration_rate..)
p3.setMigrationRates(p2, ..migration_rate..)
p2.tag = 1;
}
if(p2.tag == 1 & abund[1, gen - 1] <= ..abund_threshold..) {
p1.setMigrationRates(p2, 0)
p3.setMigrationRates(p2, 0)
p2.tag = 0;
}
if(p3.tag == 0 & abund[2, gen - 1] > ..abund_threshold..) {
p1.setMigrationRates(p3, ..migration_rate..)
p2.setMigrationRates(p3, ..migration_rate..)
p3.tag = 1;
}
if(p3.tag == 1 & abund[2, gen - 1] <= ..abund_threshold..) {
p1.setMigrationRates(p3, 0)
p2.setMigrationRates(p3, 0)
p3.tag = 0;
}
## output Fst
slimr_output(c(calculateFST(p1, p2), calculateFST(p1, p3), calculateFST(p2, p3)), "fsts");
}),
slim_block(1000, late(), {
sim.simulationFinished()
})
) -> pop_sim_fstNow let’s run that for our defaults to see how it works.
fst_script <- slim_script_render(pop_sim_fst)## Warning: Warning: A templated variable was not specified in the template and has been replaced by its default value.fst_res <- slim_run(fst_script, throw_error = TRUE)##
##
## Simulation finished with exit status: 0##
## Success!So now we can have a look at how the Fst data was returned, which is the default way slimr returns data that is a vector:
fst_res$output_data## # A tibble: 1,000 x 5
## generation name expression type data
## <int> <chr> <chr> <chr> <chr>
## 1 1 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.0~ 0.0440576 0~
## 2 2 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.0~ 0.0479792 0~
## 3 3 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.0~ 0.0551762 0~
## 4 4 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.0~ 0.0594098 0~
## 5 5 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.0~ 0.0689084 0~
## 6 6 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.0~ 0.0978984 0~
## 7 7 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.1~ 0.105959 0.~
## 8 8 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.1~ 0.113839 0.~
## 9 9 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.1~ 0.113227 0.~
## 10 10 fsts c(calculateFST(p1, p2), calc~ "float [0:2] 0.1~ 0.114101 0.~
## # ... with 990 more rows
typeof(fst_res$output_data$data)## [1] "character"So it is currently in the form of a character vector separated by spaces. This is pretty straightforward to extract, so let’s do it (Note: slimr has the ability to automatically convert this output to more user friendly form, but here we will do it the hard way to show how you can flexibly extract different kinds of data, even if a built-in slimr function doesn’t exist to handle the type of output you want).
extract_fst <- function(fst_res) {
fst_dat <- fst_res$output_data %>%
dplyr::select(generation, data) %>%
dplyr::mutate(data = str_split(data, " ")) %>%
tidyr::unnest_longer(data, values_to = "fst", indices_to = "index") %>%
dplyr::mutate(fst = as.numeric(fst), `subpop\npair` = c("BL-BR", "BL-TR", "BR-TR")[index])
fst_dat
}
fst_dat <- extract_fst(fst_res)
ggplot(fst_dat, aes(generation, fst)) +
geom_path(aes(colour = `subpop\npair`)) +
scale_y_sqrt() +
ggforce::facet_zoom(x = generation > 750 & generation < 1000,
y = fst > 0.1,
horizontal = FALSE,
zoom.size = 1)
We can see that pairwise Fst rapidly increases between all three subpopulations, reaching nearly its theoretical maximum of 1.0, and then periodically dipping slightly. Let’s look at how this equilibrium pattern corresponds with population sizes. We will look at mean population size and mean Fst.
pop_mean <- apply(pop_abunds, 2, mean) %>%
rep(length.out = n_distinct(fst_dat$generation)) * 100 ## don't forget popsize_scaling
plot_cycle <- function(fst_dat, pop_mean) {
fst_pop <- fst_dat %>%
dplyr::group_by(generation) %>%
dplyr::summarise(fst = mean(fst),
.groups = "drop") %>%
dplyr::mutate(popsize = pop_mean) %>%
dplyr::filter(generation > 750) ## only take data after equilibrium reached
ggplot(fst_pop, aes(popsize, fst)) +
geom_path(aes(colour = generation),
arrow = arrow(length = unit(0.15, "cm"), type = "closed"),
alpha = 0.5) +
annotate("point", x = fst_pop$popsize[1], fst_pop$fst[1], colour = "red", size = 2) +
scale_x_log10() +
scale_colour_viridis_c() +
theme_minimal()
}
plot_cycle(fst_dat, pop_mean)
The plot above shows how Fst and population size are changing through time together for 250 generations (the red dot is the starting point at generation 750). It appears that as population size begins increasing the Fst begins to drop slightly. This drop in Fst continues even as population size begins to decrease again (a lag effect?), then Fst suddenly increases back to nearly 1.0 as population size reaches a low. I suspect this pattern is due to ephemeral effects of increased efficiency of selection when population size is high. As new mutations begin sweeping to fixation there is an intermediate period where they are at intermediate frequency – this has the tendency to lead to reductions in within-population heterozygosity (relative to total heterozygosity), and leads to the subpopulation becoming slightly more similar due to chance collisions between genes being selected in both populations.
Due to overplotting however, some of the details of the plot may be difficult to see. Instead let’s use gganimate to create an animated version that might show the dynamics a bit better.
animate_cycle <- function(fst_dat, pop_mean) {
fst_pop <- fst_dat %>%
dplyr::group_by(generation) %>%
dplyr::summarise(fst = mean(fst), .groups = "drop") %>%
dplyr::mutate(popsize = pop_mean) %>%
dplyr::filter(generation > 250)
anim <- ggplot(fst_pop, aes(popsize, fst)) +
geom_point() +
scale_x_log10() +
transition_reveal(generation) +
shadow_wake(0.02, falloff = "sine-in") +
ggtitle("Generation: {frame_along}") +
theme_minimal()
animate(anim, nframes = n_distinct(fst_pop$generation) * 5, fps = 30,
start_pause = 10, end_pause = 10, renderer = gifski_renderer())
}
animate_cycle(fst_dat, pop_mean)
The above animation shows 750 generations of the dynamics, with the black dot showing the position of the population in terms of the mean population size, and the mean pairwise Fst, at different time points. We can clearly see the cycles occurring and its direction.
In any event, this is obviously a pretty unrealistic scenario, and it bears little resemblance to what is happening in our actual data. So, let us try a few other sets of parameters, and see if we can get something closer to what we observe in our data. Now we can finally see how we specify non-default values to our templated variables using slim_script_render. Let’s start by changing the migration_rate parameter, which controls how much migration we get between subpopulations when they go over their threshold population size. Let’s set it to maximum migration (0.5) – this model is similar to our hypothesis for the pattern we see. We will also set selection to practically zero, so we can isolate drift as a mechanism. We set a parameter value in our script like this:
## set migration rate to 0.5 (maximum migration, e.g. 50% of population switches place)
## population will become panmictic over a certain popsize threshold
fst_script_2 <- slim_script_render(pop_sim_fst, template = list(migration_rate = 0.5,
selection_strength = 1e-12))## Warning: Warning: A templated variable was not specified in the template and has been replaced by its default value.## now run the sim and show plots of results
fst_res_2 <- slim_run(fst_script_2, throw_error = TRUE)##
##
## Simulation finished with exit status: 0##
## Success!fst_dat_2 <- extract_fst(fst_res_2)
ggplot(fst_dat_2, aes(generation, fst)) +
geom_path(aes(colour = `subpop\npair`)) +
scale_y_sqrt() +
ggforce::facet_zoom(x = generation > 750 & generation < 1000,
y = fst < 0.5,
horizontal = FALSE,
zoom.size = 1)
plot_cycle(fst_dat_2, pop_mean)
animate_cycle(fst_dat_2, pop_mean)
This shows a very distinct pattern relative to our last run. Now we see that Fst tends to be very low, particularly when population sizes are high (e.g. when migration is “turned on”), as expected. Fst increases in the intervening periods of low population size, which can only be explained by drift. Our Fst values are closer to the magnitude of values we see in our data, which max out at around 0.05, though they still tend to be much higher in this simulation overall.
Let’s try one more scenario before we attempt to more formally compare the simulations to our data. Let’s set our abund_threshold controlling at what population size the subpopulation begins exchanging individuals to zero. This means that migration will be on all the time. We will also increase selection.
fst_script_3 <- slim_script_render(pop_sim_fst, template = list(migration_rate = 0.5,
selection_strength = 0.2,
abund_threshold = 0))## Warning: Warning: A templated variable was not specified in the template and has been replaced by its default value.fst_res_3 <- slim_run(fst_script_3, throw_error = TRUE)##
##
## Simulation finished with exit status: 0##
## Success!fst_dat_3 <- extract_fst(fst_res_3)
ggplot(fst_dat_3, aes(generation, fst)) +
geom_path(aes(colour = `subpop\npair`)) +
scale_y_sqrt() +
ggforce::facet_zoom(x = generation > 750 & generation < 1000,
y = fst < 0.2,
horizontal = FALSE,
zoom.size = 1)
plot_cycle(fst_dat_3, pop_mean)
animate_cycle(fst_dat_3, pop_mean)